library(plyr)
library(DataCombine)
library(estimatr)
library(fixest)
library(plm)
library(did)
library(broom)
library(texreg)
library(usmap)
library(gsynth)
library(cspp)
library(ggplot2)

setwd("./data")



############Right to Work################

data <- readRDS("sdi_data.RDS")

rtw <- read.csv("right_to_work.csv")

unionmem <- read.csv("state_unionmembership_1983_2023.csv")
unionmem <- unionmem[unionmem$sector=="Total", c("state", "year", "members_rate", "coverage_rate")]
unionmem$year <- as.numeric(unionmem$year)
unionmem$members_rate <- as.numeric(unionmem$members_rate)*100
unionmem$coverage_rate <- as.numeric(unionmem$coverage_rate)*100



cspp <- get_cspp_data(vars=c("ogdpman", "gsptotal", "poptotal"), years=seq(2000, 2018))
cspp$ogdpman_share <- 100*(cspp$ogdpman/cspp$gsptotal)
cspp[,c("state_fips", "stateno", "state_icpsr")] <- NULL

fips <- read.csv("state_names_fips.csv")
fips <- fips[fips$st!="DC",]

data <- join(fips, data)
data <- join(data, rtw)
data <- join(data, cspp)
data <- join(data, unionmem)

data$rtw <- 0
data$rtw[!is.na(data$rtw_year) & data$year>=data$rtw_year] <- 1

data$rtw[data$st=="MO"] <- 0

data$rtw_plot <- NA




alexhf <- read.csv("statecapture_chp6_troika_antilabor.csv")
alexhf <- join(alexhf, data[data$year==2018, c("year","st","rtw")])
alexhf$future_rtw <- alexhf$rtw
alexhf$future_rtw[alexhf$rtw==1 & (alexhf$st %in% c("WV","KY","IN","MI","WI"))==F] <- NA

pdf("conservative_troika.pdf", w=4, h=4)
ggplot(alexhf, aes(x=Troika.Index, y=future_rtw)) +
  geom_point() +
  geom_smooth(method="lm", color="black") +
  ylab("Later Adopts RTW") +
  xlab("Conservative Troika Score") +
  theme_classic()
dev.off()

data <- join(data, read.csv("statecapture_chp6_troika_antilabor.csv"))

robustness_check_troika <- lm_robust(democracy_mcmc ~ rtw + unified_r + Troika.Index + pct_black + pct_latino + ogdpman_share + log(poptotal) + factor(state) + factor(year), data, clusters = state)



for(i in levels(factor(data$state))){
  
  if(length(data$state[data$state==i & data$rtw==1])==19){
    
    data$rtw_plot[data$state==i] <- "Always RTW"
    
  }
  
  if(length(data$state[data$state==i & data$rtw==0])==19){
    
    data$rtw_plot[data$state==i] <- "Never RTW"
    
  }
  
}

data$rtw_plot[is.na(data$rtw_plot)] <- as.character(data$rtw[is.na(data$rtw_plot)])
data$st_plot[data$rtw_plot %in% c("0","1")] <- data$st[data$rtw_plot %in% c("0","1")]
#data$st_plot[is.na(data$st_plot)] <- "Other"

data$rtw_plot[data$rtw_plot=="0"] <- "Before RTW"
data$rtw_plot[data$rtw_plot=="1"] <- "After RTW"


#Map

plot_means <- data[data$year==2018,]
plot_means$region <- tolower(plot_means$state)

plot_means$rtw_year[plot_means$st=="MO"] <- NA

plot_means$full <- plot_means$state

data(statepop)

statepop <- join(statepop, plot_means)
statepop$region <- NULL

# Get centroids
centroid_labels <- usmapdata::centroid_labels("states")

# Join data to centroids
data_labels <- merge(centroid_labels, statepop, by = "fips")


pdf("rtw_map.pdf")
plot_usmap(data = statepop, values = "rtw", color = "black", labels = FALSE) +
  guides(fill = "none") +
  scale_fill_gradient(low = "white", high = "grey", name="") +
  geom_text(data = data_labels, ggplot2::aes(x = x, y = y, label = rtw_year), color = "black", size=2.75)
dev.off()





###Granger test

data <- slide(data, Var="democracy_mcmc", TimeVar="year", GroupVar="state", NewVar="democracy_mcmc_lag",
              slideBy=-1)
m1 <- lm_robust(rtw ~ democracy_mcmc_lag + factor(state) + factor(year), data, clusters = state)

m2 <- lm_robust(rtw ~ democracy_mcmc_lag + unified_r + factor(state) + factor(year), data, clusters = state)

data$democracy_mcmc_lag <- NULL
data <- slide(data, Var="democracy_mcmc", TimeVar="year", GroupVar="state", NewVar="democracy_mcmc_lag",
              slideBy=-2)
m3 <- lm_robust(rtw ~ democracy_mcmc_lag + factor(state) + factor(year), data, clusters = state)

m4 <- lm_robust(rtw ~ democracy_mcmc_lag + unified_r + factor(state) + factor(year), data, clusters = state)

data$democracy_mcmc_lag <- NULL
data <- slide(data, Var="democracy_mcmc", TimeVar="year", GroupVar="state", NewVar="democracy_mcmc_lag",
              slideBy=-3)
m5 <- lm_robust(rtw ~ democracy_mcmc_lag + factor(state) + factor(year), data, clusters = state)

m6 <- lm_robust(rtw ~ democracy_mcmc_lag + unified_r + factor(state) + factor(year), data, clusters = state)

data$democracy_mcmc_lag <- NULL
data <- slide(data, Var="democracy_mcmc", TimeVar="year", GroupVar="state", NewVar="democracy_mcmc_lag",
              slideBy=-4)
m7 <- lm_robust(rtw ~ democracy_mcmc_lag + factor(state) + factor(year), data, clusters = state)

m8 <- lm_robust(rtw ~ democracy_mcmc_lag + unified_r + factor(state) + factor(year), data, clusters = state)


texreg(list(m1, m2, m3, m4, m5, m6, m7, m8), include.ci = F, stars = numeric(0), 
       digits=4, omit.coef = c("factor|(Intercept)"), booktabs=T, 
       include.rsquared=F, include.adjrs=F,
       file="granger_results.tex")


####TWFE



m1 <- lm_robust(democracy_mcmc ~ rtw + factor(state) + factor(year), data, clusters = state)

m2 <- lm_robust(democracy_mcmc ~ rtw + unified_r + factor(state) + factor(year), data, clusters = state)

m3 <- lm_robust(democracy_mcmc ~ rtw + pct_black + pct_latino + ogdpman_share + log(poptotal) + factor(state) + factor(year), data, clusters = state)

m4 <- lm_robust(democracy_mcmc ~ rtw + unified_r + pct_black + pct_latino + ogdpman_share + log(poptotal) + factor(state) + factor(year), data, clusters = state)


texreg(list(m1, m2, m3, m4), include.ci = F, stars = numeric(0), 
       digits=4, omit.coef = c("factor|(Intercept)"), booktabs=T, 
       include.rsquared=F, include.adjrs=F,
       file="twfe_results.tex")


#union membership

m1 <- lm_robust(democracy_mcmc ~ members_rate + factor(state) + factor(year), data, clusters = state)

m2 <- lm_robust(democracy_mcmc ~ members_rate + unified_r + factor(state) + factor(year), data, clusters = state)

m3 <- lm_robust(democracy_mcmc ~ members_rate + pct_black + pct_latino + ogdpman_share + log(poptotal) + factor(state) + factor(year), data, clusters = state)

m4 <- lm_robust(democracy_mcmc ~ members_rate + unified_r + pct_black + pct_latino + ogdpman_share + log(poptotal) + factor(state) + factor(year), data, clusters = state)

texreg(list(m1, m2, m3, m4), include.ci = F, stars = numeric(0), 
       digits=4, omit.coef = c("factor|(Intercept)"), booktabs=T, 
       include.rsquared=F, include.adjrs=F,
       file="twfe_results_membership.tex")




##Descriptive Plots

data$rtw_plot <- factor(data$rtw_plot, levels = c("Before RTW", "After RTW", "Always RTW", "Never RTW"))


pdf("rtw_descriptive.pdf", h=6, w=8)
ggplot(data, aes(x=year, y=democracy_mcmc, color=rtw_plot, linewidth=rtw_plot)) +
  geom_line(aes(group=st)) +
  geom_text(aes(label=st_plot), color="black") +
  #geom_point(data=data[data$rtw_plot %in% c("Before RTW", "After RTW"),], aes(x=year, y=democracy_mcmc, shape=st), size=2, color="black") +
  #geom_smooth(aes(group=factor(partycontrol))) +
  scale_color_manual(values=c("dodger blue","red","grey60","grey90"), name="") +
  scale_shape_manual(values=c(1,2,3,4,5,6), guide=F) +
  scale_linewidth_manual(values=c(0.75, 0.75, 0.5, 0.5), name="", guide=F) +
  xlab("Year") +
  ylab("State Democracy Score") +
  scale_x_continuous(breaks=c(2000,2005,2010,2015)) +
  theme_classic()
dev.off()


data$rtw_plot_bw <- data$rtw_plot
data$rtw_plot_bw[data$rtw_plot]

pdf("rtw_descriptive_bw.pdf", h=6, w=8)
ggplot(data[data$rtw_plot!="Always RTW" & data$rtw_plot!="Never RTW",], 
       aes(x=year, y=democracy_mcmc, color=rtw_plot, shape=rtw_plot, linewidth=rtw_plot)) +
  geom_line(aes(group=st)) +
  geom_text(aes(label=st_plot), color="black") +
  #geom_point() +
  #geom_point(data=data[data$rtw_plot %in% c("Before RTW", "After RTW"),], aes(x=year, y=democracy_mcmc, shape=st), size=2, color="black") +
  #geom_smooth(aes(group=factor(partycontrol))) +
  scale_color_manual(values=c("grey","black"), name="") +
  #scale_shape_manual(values=c(1,2,3,4)) +
  scale_linewidth_manual(values=c(0.75, 0.75, 0.5, 0.5), name="", guide=F) +
  xlab("Year") +
  ylab("State Democracy Score") +
  scale_x_continuous(breaks=c(2000,2005,2010,2015)) +
  theme_classic()
dev.off()




periods <- data.frame(seq(2000, 2018, by=1), c(1:19))
names(periods) <- c("year", "period")

data$period <- NULL
data <- join(data, periods, match="all")

data$first_treat_period <- NA
data$invariant_r <- 0


for(i in levels(factor(data$st))){
  
  data$first_treat_period[data$st==i] <- min(data$period[data$st==i & data$rtw==1], na.rm=T)
  
  if(data$unified_r[data$st==i & data$year==2011]==1 & !is.na(data$unified_r[data$st==i & data$year==2011])){
    
    data$invariant_r[data$st==i] <- 1
    
  }
}

data$first_treat_period[is.infinite(data$first_treat_period)|is.na(data$first_treat_period)] <- 0

data$id <- as.numeric(as.factor(data$st))


att <- att_gt(yname = "democracy_mcmc",
              tname = "period",
              idname = "id",
              gname = "first_treat_period",
              #xformla = ~X,
              allow_unbalanced_panel = F,
              data = data
)


att_simple <- aggte(att, type = "simple")
att_group <- aggte(att, type = "group")
att_dynamic <- aggte(att, type = "dynamic")
att_calendar <- aggte(att, type = "calendar")

outlist <- list()

for(i in ls(pattern="att_")){
  
  temp <- get(i)
  temp_out <- data.frame(temp$overall.att, temp$overall.se, i)
  names(temp_out) <- c("estimate","se","model")
  
  outlist[[i]] <- temp_out
  
}

plotdata <- do.call(rbind, outlist)
plotdata$ATT <- c("Calendar", "Dynamic", "Group", "Group-Time")

plotdata$cilow <- plotdata$estimate - 1.96*plotdata$se
plotdata$cihigh <- plotdata$estimate + 1.96*plotdata$se


pdf("rtw_callawaysantanna.pdf",h=6,w=6)
ggplot(plotdata, aes(x=ATT, y=estimate)) +
  geom_point() +
  geom_errorbar(aes(ymin=cilow, ymax=cihigh), width=0) +
  ylab("Effect on State Democracy Score\n(in Standard Deviations)") +
  xlab("ATT Aggregation") +
  geom_hline(yintercept=0, linetype=2) +
  theme_bw()
dev.off()


######Event Study

did_data <- data

did_data$treat_period_num <- did_data$period - did_data$first_treat_period

did_data$ever_treated[did_data$first_treat_period!=0] <- 1
did_data$ever_treated[did_data$first_treat_period==0] <- 0

did_data$always_treated <- 0
did_data$always_treated[did_data$rtw_year<2000] <- 1

did_data$treat_period <- did_data$treat_period_num


did_data$treat_period <- did_data$treat_period*did_data$ever_treated

did_data$treat_period <- factor(did_data$treat_period)

did_data$treat_period <- relevel(factor(did_data$treat_period), ref = 17)


temp <- estimatr::tidy(lm_robust(democracy_mcmc ~ treat_period + factor(state) + factor(period),
                                 se_type = "stata",
                                 clusters = state,
                                 did_data[!is.na(did_data$first_treat_period) & did_data$always_treated==0,]))

temp <- temp[grepl("treat_period", temp$term),]
temp_baseline <- temp[1,]
temp_baseline$term <- "treat_period-1"
temp_baseline$estimate <- 0
temp_baseline$std.error <- NA
temp <- rbind(temp, temp_baseline)
temp$model <- "Full RTW Effect"


temp_cov <- estimatr::tidy(lm_robust(democracy_mcmc ~ treat_period + unified_r + factor(state) + factor(period),
                                 se_type = "stata",
                                 clusters = state,
                                 did_data[!is.na(did_data$first_treat_period) & did_data$always_treated==0,]))

temp_cov <- temp_cov[grepl("treat_period", temp_cov$term),]
temp_cov_baseline <- temp_cov[1,]
temp_cov_baseline$term <- "treat_period-1"
temp_cov_baseline$estimate <- 0
temp_cov_baseline$std.error <- NA
temp_cov <- rbind(temp_cov, temp_cov_baseline)
temp_cov$model <- "RTW Effect Net of GOP Control"

plotdata <- rbind(temp, temp_cov)

plotdata$term_clean <- as.numeric(as.character(gsub("treat_period", "", plotdata$term)))
plotdata$cilow <- plotdata$estimate - (1.96*plotdata$std.error)
plotdata$cihigh <- plotdata$estimate + (1.96*plotdata$std.error)



pdf("eventstudy_rtw.pdf", h=6, w=12)
ggplot(plotdata[plotdata$term_clean>=-8 & plotdata$term_clean<=8,], aes(x=term_clean, y=estimate)) +
  geom_errorbar(aes(ymin=cilow, ymax=cihigh), width=0.5) +
  #geom_errorbar(aes(ymin=cilow, ymax=cihigh), width=0, color="grey") +
  #geom_point() +
  #geom_line() +
  geom_point() +
  #geom_segment(data=lines_pre[lines_post$model=="Controls",], aes(x = -16, y = estimate, xend = 0, yend = estimate)) +
  #geom_segment(data=lines_post[lines_post$model=="Controls",], aes(x = 0, y = estimate, xend = 16, yend = estimate)) +
  facet_wrap(~model) +
  scale_x_continuous(breaks=c(-12,-8,-4,0,4,8,12)) +
  #scale_y_continuous(breaks=c(-0.1,-0.05,0,0.05,0.1), limits=c(-.1,.1)) +
  geom_hline(yintercept=0, linetype=2) +
  geom_vline(xintercept=0, linetype=2) +
  xlab("Years Until RTW") +
  ylab("Effect of RTW on State Electoral Democracy") +
  theme_bw()
dev.off()







############Leave-One-Out Analysis##############



periods <- data.frame(seq(2000, 2018, by=1), c(1:19))
names(periods) <- c("year", "period")

data$period <- NULL
data <- join(data, periods, match="all")

data$first_treat_period <- NA
data$invariant_r <- 0


for(i in levels(factor(data$st))){
  
  cat(i, "\n")
  
  data$first_treat_period[data$st==i] <- min(data$period[data$st==i & data$rtw==1], na.rm=T)
  
  if(data$unified_r[data$st==i & data$year==2011]==1 & !is.na(data$unified_r[data$st==i & data$year==2011])){
    
    data$invariant_r[data$st==i] <- 1
    
  }
}

data$first_treat_period[is.infinite(data$first_treat_period)|is.na(data$first_treat_period)] <- 0

data$id <- as.numeric(as.factor(data$st))

outlist2 <- list()

for(j in levels(factor(data$st))){
  
  att <- att_gt(yname = "democracy_mcmc",
                tname = "period",
                idname = "id",
                gname = "first_treat_period",
                #xformla = ~X,
                allow_unbalanced_panel = F,
                data = data[data$st!=j,]
  )
  
  
  att_simple <- aggte(att, type = "simple")
  att_group <- aggte(att, type = "group")
  att_dynamic <- aggte(att, type = "dynamic")
  att_calendar <- aggte(att, type = "calendar")
  
  outlist <- list()
  
  for(i in ls(pattern="att_")){
    
    temp <- get(i)
    temp_out <- data.frame(temp$overall.att, temp$overall.se, i)
    names(temp_out) <- c("estimate","se","model")
    
    outlist[[i]] <- temp_out
    
  }
  
  plotdata <- do.call(rbind, outlist)
  plotdata$ATT <- c("Calendar", "Dynamic", "Group", "Group-Time")
  
  plotdata$cilow <- plotdata$estimate - 1.96*plotdata$se
  plotdata$cihigh <- plotdata$estimate + 1.96*plotdata$se
  
  plotdata$leftout <- j
  
  outlist2[[j]] <- plotdata
  
}

plotdata <- do.call(rbind, outlist2)


plotdata$leftout <- factor(plotdata$leftout, levels=rev(levels(factor(plotdata$leftout))))

pdf("rtw_callawaysantanna_leaveoneout.pdf",h=7,w=7)
ggplot(plotdata, aes(x=leftout, y=estimate)) +
  geom_point() +
  geom_errorbar(aes(ymin=cilow, ymax=cihigh), width=0) +
  facet_wrap(~ATT, ncol=4) +
  coord_flip() +
  ylab("Effect on State Democracy Score\n(in Standard Deviations)") +
  xlab("State Left Out") +
  geom_hline(yintercept=0, linetype=2) +
  theme_bw()
dev.off()






#Leaving out IN and MI

att <- att_gt(yname = "democracy_mcmc",
              tname = "period",
              idname = "id",
              gname = "first_treat_period",
              #xformla = ~X,
              allow_unbalanced_panel = F,
              data = data[(data$st %in% c("IN", "MI"))==F,]
)


att_simple <- aggte(att, type = "simple")
att_group <- aggte(att, type = "group")
att_dynamic <- aggte(att, type = "dynamic")
att_calendar <- aggte(att, type = "calendar")

outlist <- list()

for(i in ls(pattern="att_")){
  
  temp <- get(i)
  temp_out <- data.frame(temp$overall.att, temp$overall.se, i)
  names(temp_out) <- c("estimate","se","model")
  
  outlist[[i]] <- temp_out
  
}

plotdata <- do.call(rbind, outlist)
plotdata$ATT <- c("Calendar", "Dynamic", "Group", "Group-Time")

plotdata$cilow <- plotdata$estimate - 1.96*plotdata$se
plotdata$cihigh <- plotdata$estimate + 1.96*plotdata$se










































